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Abstract — In this paper, we improve iterative greedy search 
algorithms in which atoms are selected serially over iterations, 
i.e., one-by-one over iterations. For serial atom selection, we 
devise two new schemes to select an atom from a set of potential 
atoms in each iteration. The two new schemes lead to two new 
algorithms. For both the algorithms, in each iteration, the set 
of potential atoms is found using a standard matched filter. In 
case of the first scheme, we propose an orthogonal projection 
strategy that selects an atom from the set of potential atoms. 
Then, for the second scheme, we propose a look ahead strategy 
such that the selection of an atom in the current iteration has 
an effect on the future iterations. The use of look ahead strategy 
requires a higher computational resource. To achieve a trade- 
off between performance and complexity, we use the two new 
schemes in cascade and develop a third new algorithm. Through 
experimental evaluations, we compare the proposed algorithms 
with existing greedy search and convex relaxation algorithms. 

Index Terms — Sparse signal estimation, compressive sampling, 
orthogonal matching pursuit, orthogonal least squares, subspace 
pursuit. 

I. Introduction 

A linear underdetermined system model based sparse sig- 
nal estimation problem has attracted much attention in the 
current literature. The sparse signal estimation problem has 
many applications, for example, in linear regression [1], 
communication [2] -[4], multimedia [5] -[8], and recently in 
compressive sampling (CS) [9]-[10]. The algorithms proposed 
for solving such a problem may be categorized into three broad 
classes: convex relaxation, [9]-[18], Bayesian inference [20]- 
[24], and iterative greedy search [25]-[41]. The iterative greedy 
search (IGS) algorithms are of rising interest for solving large 
dimensional sparse signal estimation problems due to their 
algorithmic simplicity and lower complexity. 

Generally an IGS algorithm constructs a support set of 
the underlying sparse signal vector through iterations. In the 
literature, the support set is referred to as the set of indices cor- 
responding to coordinates of non-zero elements of the sparse 
signal vector. Also the columns of the measurement matrix (or 
dictionary matrix) are referred to as atoms. Naturally, the sup- 
port set's elements are the indices of atoms associated with the 
non-zero elements of the sparse vector. In the underdetermined 
setup, an atom-selection process for constructing the support 
set is a critical task. Depending on atom selection approaches, 
the IGS algorithms may be categorized into two classes: 
serial (or sequential) and parallel (or simultaneous). Prominent 
examples of the serial atom selection based IGS algorithms are 
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orthogonal matching pursuit (OMP) [26]-[27] and orthogonal 
least squares (OLS) [28]-[30]. For both the algorithms, a single 
atom is selected from a set of potential atoms in each iteration 
and then the atom's index is inducted into the support set under 
progressive construction. Therefore, the cardinality of the 
support set under construction is increased one-by-one through 
iterations serialljQ. At the end of all iterations, the algorithms 
provide a support set of full cardinality. On the other hand, 
for the parallel atom selection based IGS algorithms, such as 
subspace pursuit (SP) [35] and CoSaMP [37], a support set 
with a fixed (pre-determined) cardinality is estimated in each 
iteratior|§. A fixed number of atoms is selected simultaneously 
in each iteration to form the support set of full cardinality and 
the estimate of the underlying support set is refined through 
iterations. In practice, IGS algorithms provide different support 
set estimation performances (in turn different sparse signal 
estimation performances). For selecting (or detecting) an atom 
in a computationally efficient manner, the use of a matched 
filter remains prevalent in existing IGS algorithms. 

In this paper, we develop new IGS algorithms to provide a 
better estimate of the underlying support set at the expense 
of computational complexity. Our interest is to deal with 
a highly underdetermined system that may arise in many 
applications, for example in a highly under-sampled CS. For 
such a underdetermined system, the reliability of the matched 
filter degrades. At the expense of complexity, we endeavor for 
improving atom selection performance beyond the scope of 
a matched filter. In the new algorithms, a trade-off between 
complexity and performance can be achieved by adjusting user 
defined tunable parameters. 

We develop the new algorithms using the architectures of 
OMP and OLS algorithms where atom selection is performed 
serially over iterations. For the new algorithms, in each itera- 
tion, a set of potential atoms is chosen and then a single atom 
is finally selected from the set of potential atoms. We propose 
new atom selection strategies and incorporate them in the OMP 
and OLS architectures. Strategically, the new algorithms use a 
two stage mechanism: (a) first, the standard matched filter is 
used to choose a set of potential atoms from all the competing 
atoms, and (b) second, an atom is finally selected from the set 
of potential atoms by using the new selection strategies. In the 
first stage, the cardinality of the set of potential atoms is pre- 

1 Whilst both the OMP and OLS algorithms are structurally similar, there 
exists a subtle difference between them in the way an atom is selected in each 
iteration. A clear exposition on the difference between these two algorithms 
can be found in [31]. The reader is encouraged to see Appendix lAl 

2 We briefly describe serial atom selection based OMP and OLS algorithms 
in Appendix fA] and parallel atom selection based SP in Appendix [E] 
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determined. The use of matched filter in the first stage helps 
to reduce search space for the new selection strategies invoked 
in the second stage. 

For the standard OMP, in each iteration, the matched filter 
is used to select an atom from all the competing atoms. In 
a highly underdetermined setup, atoms are more correlated 
and hence the atom selection performance using matched filter 
degrades. This performance degradation is due to fact that an 
underlying atom is subjected to a considerable interference 
from other closely correlated atoms. To improve atom selec- 
tion, we propose to use the same matched filter for choosing 
a set of potential atoms and then select an atom from the set 
using an orthogonal projection strategy. The new algorithm is 
referred to as the projection based OMP (POMP) where the 
use of orthogonal projection is motivated by its existing use 
in parallel atom selection based IGS algorithms (such as in 
SP and CoSaMP). The orthogonal projection invokes signal 
estimation in least-squares (LS) sense and hence its use adds 
more strength on the top of matched filter. 

Next, for the standard OLS, an atom is selected in each 
iteration such that the selection leads to the minimum norm of 
a fitting residual (LS residual) for that iteration. The selection 
of an atom in the current iteration does not depend on the 
final performance of the OLS algorithm when all the iterations 
would have been finished. Therefore, the selection of a new 
atom in the current iteration is blind to its effect on the future 
iterations in the sense of minimizing the final fitting residual 
norm. To overcome the shortcoming of blindness in OLS, we 
propose to use a look ahead selection strategy [39]. In the 
current iteration, an atom is selected from the set of potential 
atoms by evaluating its effect on the final performance measure 
in the sense of minimizing the norm of fitting residual at the 
end of all future iterations. The new algorithm is referred to 
as look ahead OLS (LAOLS). Use of a look ahead strategy is 
shown to provide a better performance, but at the expense of 
a higher computation. 

Further, we combine the orthogonal projection-based and 
look ahead atom selection strategies in cascade for developing 
an IGS algorithm which provides a trade-off between compu- 
tational complexity and reconstruction performance. 

Considering compressive sampling (CS) as the application, 
we experimentally evaluate the new algorithms and compare 
their performances vis-a-vis several existing algorithms, such 
as OMP, OLS, SP and convex relaxation algorithms. For 
experiments, we endeavor for CS reconstruction of Gaussian 
and binary sparse signals at varying under-sampling rates and 
measurement noise levels. 

Notations: Let A E R MxN , x e R N , and X C 
{1, 2, . . . , N}. We denote \X\ and X as the cardinality and 
complement of the set X, respectively. The matrix Aj G 
R.Mx|z| consists of the columns of A indexed by i € X, and 
Xi € R' 1 ' is composed of the components of x indexed by 
i e I. Also we denote (.)' and {.)' as transpose and pseudo- 
inverse, respectively. 



A. Preliminaries of Compressive Sampling 

Let us state the standard CS problem where we acquire a 
If -sparse signal x £ M. N via the linear measurements 

y = Ax + w, (1) 

where A g IR MxAr is a matrix representing the sampling 
system, y £ M. represents a vector of measurements and 
w G K. is additive noise representing measurement errors. A 
If -sparse signal vector consists of at most K non-zero scalar 
components. For the setup K < M < N (underdetermined 
system of linear equations), the task is to reconstruct x from 
y as x. With no a-priori statistical knowledge of x and w, 
the objective in CS is to strive for a reduced number of 
measurements (M) as well as achieving a good reconstruction 
quality. Note that, in practice, we may wish to acquire a 
signal x that is sparse in a known orthonormal basis and the 
concerned problem can be recast as fl}. 

For a sparse signal vector x = [x\, X2, ■ ■ ■ , Xn] , the 
support set I C {1,2,..., N} is defined as I = {i : xi ^ 0}. 
For a X-sparse vector x € Mr, \X\ — ||x||o < K. In this 
paper, we assume that \X\ = K. Denoting the z'th column 
(atom) of the measurement matrix A as a s , note from (Q~|) that 

y = ^]s ! a i +w = Aixi + w, (2) 

iei 

where xj € 1^ and Ax € K M x K ■ From y, if the underlying 
support set X can be estimated, then we can estimate the 
non-zero values of x (i.e. xj) using a standard least squares 
(LS) solution (as K < M, we can use the pseudo-inverse). 
Therefore, a better estimate of the support set leads to a better 
reconstruction performance. However, if the matrix Ax is not 
full column rank (or ill-conditioned) then the LS solution 
will be erroneous. This may happen in case of the matrices 
where the correlation between the columns are considerably 
high, for example, in a highly under-sampled CS where 
M <C N. Recently, theoretical properties of such matrices 
are investigated in [42], [43]. 

Next we consider the issue of constructing the sensing 
matrix A G M. MxN . While it is possible to obtain deterministic 
constructions of A = {ay} holding a specific structure, at 
present the most efficient designs (i.e., those requiring min- 
imum number of rows) rely on random matrix constructions 
where the ay's are assumed realizations of independent and 
identically distributed (i.i.d.) random variables. A standard 
method is to draw ay's independently from a Gaussian source 
(i.e., ay ~ N (0, jr)) and then to scale the columns of A to 
unit-norm. Note that once A is constructed, it remains fixed 
and made known to a CS reconstruction algorithm. In the 
following sections, we develop the new IGS algorithms. 

II. Projection Based Orthogonal Matching 
Pursuit 

Using the algorithmic structure of OMP, we develop projec- 
tion based OMP (POMP) where the selection of an atom in the 
current iteration is carried out using an orthogonal projection 
method. 
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For developing the POMP algorithm, we first need to 
describe the orthogonal projection based selection method in 
Algorithm [T] In Algorithm Q] the inputs are A, y, previous 



Algorithm 1 : Projection-Based Atom-index Selection 
Input: 
1: A, y; 

2: Previous support set: I; 
3: Set of potential atoms' indices: 
Assumption: 

i: i (p) ni = 0, x e R N ; 

Execution: 

2: xj (u) Aj ( } y; (Orthogonal projection) 

3: Xj ( ^ <— 0; (Zeroing the coordinates indexed by 2~( p )) 
4: i 4— index of the highest amplitude of x 
Output: 

1: The index i. (Note: i € I( p )) 



support set X and the set of potential atoms' indices 
We refer to the previous support set I as the constructed 
support set in the previous iteration. The T/ p \ is the set of 
potential atoms chosen using the standard matched filter in the 
current iteration. In Algorithm Q] a dummy set is formed as 
X(„) = IU1m. Then the measurement vector y is projected 
on Ax, , followed by constructing an estimate of x where 
} = 0. The index of the most promising atom is then 
found using the standard criterion of finding the coordinate 
associated with the highest amplitude of x. Such an orthogonal 
projection based atom selection mechanism is used in parallel 
IGS algorithms, for example in SP [35] and CoSaMP [37], 
albeit differently. To compare, the SP algorithm is described 
briefly in Appendix IB1 

Using Algorithm [T] let us define the following function. 

Function 1: (Projection-based atom-index selection) Let 
y € R M , A e R MxN , previous support set 1 and the set 
of potential atoms' indices 2~( p ) are given. Suppose i denotes 
the index of the new atom selected in the current iteration. 
Then i is the output of the following algorithmic function 

i = proj_atom_index (A, y , I, Ir p ) ) (3) 

where the above function executes the Algorithm Q] 

Using Function [T] the main steps of the POMP algorithm 
are summarized in Algorithm |2] where the support set is 
constructed serially. In the POMP algorithm, L is a user 
defined positive integer parameter that decides how many 
atoms need to be checked as potential atoms in each iteration. 
The algorithm starts with an initial empty support set Iq = 
and an initial residual ro = y. At the fc'th iteration stage, 
it uses the standard matched filter to choose the set of L 
potential atoms' indices (in step 0, orthogonal projection 
method to select a new atom's index from the set (in stepHJ), 
inducts the selected index into the intermediate support set 
under construction (in step [5]), solves a LS problem with 
the intermediate support set and produces a new residual by 
subtracting the LS fit (in step |6j. Given the sparsity level 



Algorithm 2 : POMP for CS Recovery 
Input: 

1: A, y, K; 

2: L (L potential atoms and 1 < L < K); 

Initialization: 

1: Iteration counter k <— 0; 

2: r <- y, T <r- 0; 
Iterations: 

i: repeat 

2: k 4- k + 1; 

3: 2~( p ) <— {indices of L highest amplitudes of A*rfc_i}; 

(Assumption: chosen L indices ^ Tk-i) 
4: ifc <— proj_atom_index (A, y,2fc_i,2~( p )); 
5: 2fc 4r- 2Vi U i k ; (Note: \X k \ = k) 

6: rfe ^— y — Aj k Aj y; (Orthogonal projection) 

7: until ((||r fc || 2 > ||r fe ^|| 2 ) or (k > K)) 
8: k <— k — 1; (Previous iteration) 

Output: 

1: x e M. N , satisfying xx k — Aj fc y and x^ fc = 0. 



K, the algorithm usually executes K iterations and forms a 
support set of cardinality K at the end. Note that the support 
set cardinality is increased one-by-one by selecting a new 
atom's index in each iteration. We also have used a stopping 
criterion that ensures a decreasing trend of residual norm 
over iterations and the maximum allowable cardinality of the 
constructed support set. We assume that the sparsity level K 
is known a-priori like the cases of OMP (see Algorithm 3 of 
[33])), SP [35] and CoSaMP [37]. 

In contrast to the OMP algorithm (see Appendix |A] for a 
brief description), the new additions in the POMP algorithm 
are the use of the set of potential atoms and the orthogonal 
projection based atom selection strategy from the set. A natural 
question is why does the projection based atom selection 
method perform better than the standard matched filter. For 
understanding, let us consider an over-determined case where 
M > N and assume that the columns of A (i.e. atoms) 
are orthonormal [| Now let us consider OMP for such an 
overdetermined A and assume clean measurement condition. 
In that case, the matched filter provides exact signal com- 
ponent values and hence the use of matched filter is optimal. 
The matched filter selects all the underlying K atoms perfectly 
one-by-one over iterations. The selection order of atoms will 
follow a strictly regular path as follows: the atom associated 
with largest amplitude scalar value will be selected first and 
then the atom associated with second largest amplitude scalar 
value will be selected, and so on. That means the atom 
selection order will follow the strict path associated with the 
decreasing amplitudes of underlying non-zero scalar signal 
values. Now let us consider the underdetermined CS setup 
(or sparse signal estimation problem) where M < N and 
hence all the atoms can not be orthonormal. In that case, 
the matched filter can not provide exact signal component 

3 On this assumption, we mention that the construction of A by using i.i.d 
random variables guarantees a near-orthonormal condition by high probability. 
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values, but a mixed output and hence the use of matched 
filter is non-optimal. For any chosen atom, there exists a set 
of atoms which are highly correlated with the chosen atom; 
this set of atoms are coherent with the chosen atom. The 
mixed output of the matched filter is the outcome of the 
fact that an underlying atom is subjected to an interference 
from other atoms, mainly from the coherent atoms. The 
interference is due to non-orthogonality and hence the inner 
product operations in the matched filter do not provide reliable 
signal component values. For example, the highest amplitude 
coefficient of the output vector of the matched filter may 
not be associated with the corresponding underlying atom. 
Thus, the highest amplitude based detection strategy in the 
matched filter may not provide perfect selection of atoms over 
iterations. If M <C N, the level of interference is considerably 
high and efficiency of the matched filter degrades severely, 
leading to poor atom selection. To circumvent the problem 
partially, we propose to use the orthogonal projection method 
in POMP. In each iteration of POMP, the orthogonal projection 
helps to estimate the signal component values associated with 
the chosen potential atoms in the LS sense. In general, the 
LS approach is powerful to estimate a signal. Thus, the LS 
signal estimation followed by the highest amplitude based 
atom selection strategy has a better potential than the sole use 
of matched filter. Note that the matched filter provides a LS 
solution for a single atom, without considering contributions 
from the non-orthogonal atoms. 

Next we discuss the complexity of the POMP algorithm 
where the orthogonal projection operation is performed twice 
in each iteration. Considering that the matched filtering and 
orthogonal projection operations are the most computationally 
intensive among all the relevant computational operations, the 
POMP complexity can be seen loosely as less than twice of 
the OMP complexity. We assume that the other computational 
operations require negligible resources. To be quantitative, 
the POMP requires to perform K matched filter operations 
which is same as that of OMP. However, the POMP requires 
2K orthogonal projection operations compared to the K 
orthogonal projection operations in OMP. In the POMP, we 
used the parameter L (associated with the number of potential 
atoms) such that L < K. Our engineering perspective for 
the choice of L < K is that no more atoms, than the total 
number of atoms that will eventually be chosen (i.e., K atoms), 
need to be tested as the potential atoms in any iteration. Also, 
another consideration is that the atom selection performance 
saturates with increasing L. Through experimental evaluations 
in section IVI-CI we show that an increase in L leads to a 
better atom selection performance initially and then saturates. 
Observe that for L = 1, the performance of POMP is the same 
as the OMP. 

In the POMP, the atom selection strategy for the current 
iteration does not depend on the future iterations. In the next 
section, we propose an algorithm that sees the future. 

III. Look Ahead Orthogonal Least Squares 

In this section, we develop look ahead OLS (LAOLS) 
algorithm where the selection of an atom's index in the 



current iteration is carried out according to its future effect on 
minimizing the final residual norm. The final residual norm 
is evaluated at the end of all future iterations. For developing 
the LAOLS algorithm, we first need to realize a look ahead 
strategy in an algorithmic manner such that the future effect 
can be evaluated. Algorithm [3] shows such a look ahead 
strategy to evaluate the final residual norm. 

Algorithm 3 : Look Ahead Residual Norm 
Input: 

1: A, y, K; 

2: Previous support set 1, current atom-index i; 
Assumption: 

l: i £ 1; 
Initialization: 

1: Iteration counter k ■<— \X U i\; 

2: It^lU i, r k <- y - A Ifc A^y; 
Iterations: 

l: repeat 

2: k <- k + 1; 

3: ik index of the highest amplitude of A*rfc_i 

(Assumption: if. ^ X&_i); 
4: l k <- Z fc _i U i k ; (Note: \X k \ = k) 

5: rfc <— y — Ax k Aj y; (Orthogonal projection) 

6: until ((||r fe || 2 > ||r fc *i|| 2 ) or (k > K)) 
1: k i— k — 1; (Previous iteration) 

Output: 

1: ||r fe || 2 (where r fc = y - A Xk A' Xh y). 



The Algorithm [3] can be seen as a look ahead part of the 
OMP algorithm in the current iteration where the future atoms 
are chosen using the matched filter and orthogonal projection 
operations. Given the sparsity level K, a previous support set 
I and a current atom-index i (index of the atom chosen in 
the current iteration), the algorithm usually finds a A"-element 
support set at the end of all future iterations and returns the 
final LS residual norm. For the given previous support set I, 
the final effect of the choice of a new atom in the current 
iteration can be evaluated through using this look ahead part. 

Using Algorithm [3] let us define the following function. 

Function 2: (Look ahead residual norm) Let y € R M , A G 
RAfxiV and the sparsity level K are given. Suppose that the 
previous support set is I and the current atom-index is i. Then 
the look ahead residual norm n £ R is the output of the 
following algorithmic function 

n = look_ahead_resid_norm (A, y, K,l, i) (4) 

where the above function executes the Algorithm [3] 

Using Function 12 the main steps of the LAOLS algorithm 
are summarized in Algorithm [4] Like POMP, we fix an integer 
parameter L < K that decides how many potential atoms 
are checked in each iteration using the look ahead strategy 
of Function [2] (or Algorithm [3j. The atom with smallest look 
ahead residual norm is selected and added to the previous 
support set l k -i to form the new support set l k of cardinality 
k. Given the sparsity level K, the algorithm usually executes 
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Algorithm 4 : LAOLS for CS Recovery 
Input: 

1: A, y, K; 

2: L (L potential atoms and 1 < L < K); 

Define: 

I: J = [ji 32, -.-Jl]*, n= [ni n 2 ,...,n L ] t ; 
Initialization: 

1: Iteration counters k <— (Outer loop), I <— (Inner loop); 

2: r <- y, Tq <r- 0; 
Iterations: 
l: repeat 

2: fc «- k + 1; 

3: I( p ) «— {indices of L highest amplitudes of A*rfc_i}; 

(Assumption: chosen L indices ^ Ik-i) 
4: j^X (p) ; 

(Elements of Xi p \ are assigned to components of j) 
5: for I = 1 to L do 

6: n/ <— look_ahead_resid_norm (A, y, K,Ik~i, 

(Final residual norm if ji is selected) 
7: end for 

8: I* 4- coordinate of the lowest component of n; 
9: i k <~ jl* ; 

(Index of atom providing lowest final residual norm) 
10: X k <- U i fe ; (Note: \X k \ = jfe) 

11: r fe <- y - Ai fe A^y; (Orthogonal projection) 

12: until ((||r fe || 2 > ||r fc _i|| 2 ) or (k > K)) 
13: k 4— k — 1; (Previous iteration) 

Output: 

1: x € K w , satisfying xj fc = Aj fc y and Xj fc = 0. 



K iterations and forms a support set of cardinality K at 
the end. We mention that the LAOLS algorithm can be 
seen as a confluence of OMP and OLS algorithms. For the 
OMP connection, note the use of matched filter in the look 
ahead part to find the indices of future atoms (see step [3] of 
Algorithm [3). Also, for the OLS connection, note how the 
index of a new atom is selected in each iteration depending 
on the minimum norm of fitting residual (see steps [8] and [9] of 
Algorithm |4|. Technically, for designing a look ahead OLS, 
it is possible to use the OLS strategy of atom selection in 
Algorithm [3] instead of using the OMP strategy. However, 
we refrain from such a strategy due to a heavy increase in 
computational complexit)0. 

In contrast to the OLS algorithm (see Appendix [A] for a 
brief description), the new additions in the LAOLS algorithm 
are the use of a set of potential atoms and the look ahead atom 
selection from the set. A natural observation is that the use of 
look ahead strategy is computationally intensive. To achieve 
computational advantage, the orthogonal projection operations 
in the look ahead strategy can be performed recursively using 
block-wise matrix inversion. An exposition for performing 

4 The atom selection strategy of OLS is more computationally intensive 
than the atom selection strategy of OMP. The reader is encouraged to see 
Appendix fA] 



recursive computation is shown in Appendix ICl For each I at 
the fcth iteration, the look ahead strategy requires to perform 
(K — (k — 1)) matched filtering and orthogonal projection 
operations. Therefore, each of the total matched filtering and 
the total orthogonal projection operations in the LAOLS are 
* Ef=i Ef=i L{K - (k - 1)) = Ef =1 L{K (k 1)) = 
L(l + 2 + . . . + K) = \L(K 2 + K). Given the sparsity level 
K, the computational complexity increases linearly with the 
number of potential atoms L. A reduction in L leads to a lower 
complexity, but at an expense of CS recovery performance. 
In the next section, we develop an algorithm with a reduced 
complexity, but without much loss in performance. 

IV. Structured Orthogonal Least Squares 

In this section, we develop a structured IGS algorithm by 
combining POMP and LAOLS algorithms such that a trade- 
off between computational complexity and reconstruction per- 
formance can be established. In each iteration, we choose L 
potential atoms by the standard matched filter followed by 
reducing the cardinality of the set of potential atoms through 
an orthogonal projection based atom selection strategy. Then, 
we use the look ahead strategy to select an atom from 
the reduced set of potential atoms. The use of the reduced 
set allows for a sharp reduction in computational burden. 
However, as the cardinality reduction of the set of potential 
atoms is performed through an orthogonal projection method, 
a little performance loss can be expected. Note that an atom 
is finally selected from the set of L potential atoms through 
using the two selection schemes which are algorithmic ally in 
a cascade connection. The new algorithm is referred to as 
structured OLS (SOLS). 

To develop the SOLS algorithm, we first need to describe 
the Algorithm [5] which is a slight modification of AlgorithmQ] 
In the Algorithm |5J L' is a positive integer that decides how 



Algorithm 5 : Projection-Based Multiple Atom-indices Selec- 
tion 

Input: 

1: A, y; 

2: Previous support set: I; 

3: Set of potential atoms' indices: 

4: V (L' < \X(p)\y, 
Assumption: 

1: X(p) nl= 0, x e R N ; 
Execution: 

2: xx (u) i— At y; (Orthogonal projection) 

3: } 4— 0; (Zeroing the coordinates indexed by I( p )) 
4: X(p') indices of the L' highest amplitudes of x; 
Output: 
1: The indices' set Z( p /). 



many atoms are selected from the set of potential atoms 
Ii p ) through orthogonal projection. Considering V < \I(p)\, 
the set of L' atoms are used for further processing. For a 
given cardinality of the set I( p ), the choice of L' is again a 
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user's prerogative that depends on the trade-off between atom 
selection precision and search complexity. 

Using Algorithm [3] let us define the following function. 

Function 3: (Projection-based multiple atom-indices selec- 
tion) Let y e R M , A e R MxAr , previous support set 1 and 
the set of potential atoms' indices It p ) are given. Suppose 
I( p ') C Itp) denotes the set of indices of V atoms selected in 
the current iteration. Then Im) is the output of the following 
algorithmic function 

-^(p') = proj_multi_atom_indices (A,y,Z, (5) 

where the above function executes the Algorithm [5] 

Using Functions [2] and [3] the main steps of the SOLS 
algorithm are summarized in Algorithmic In the algorithm, we 

Algorithm 6 : SOLS for CS Recovery 
Input: 

l: A, y, K; 

2: L (L potential atoms and 1 < L < K); 

3: 7 (A subset selection parameter and < 7 < 1). 

Define: 

1: V = L — [jL\; (L' > 0; non-zero cardinality of subset) 
2: j = [ji i2,---,jL'] t , n= [m n 2 ,...,n L/ ] t ; 
Initialization: 

1: Iteration counters /c < — (Outer loop), I 4— (Inner loop); 

2: r 4-y,X 4- 0; 
Iteration: 
i: repeat 

2: k 4- k + 1; 

3: X/p\ 4— indices of the L highest amplitudes of A*r/-_i; 

(Assumption: chosen L indices ^ Ife-i) 
4: I/ p n proj_multi_atom_indices (A, y, Ife-i, 2"( p ) , £') ; 

(Reducing the number of potential atoms: \I( p i) \ = L') 

(Elements of I( P ') are assigned to components of j) 
6: for I = 1 to L' do 

7: n/ -S— look_ahcad_rcsid_norm (A, y, K,Ik-ii 
8: end for 

9: Z* 4- coordinate of the lowest component of n; 
10: ik 4- ji* ; 

(Index of the atom providing lowest final residual norm) 

li: I fe <- 2fe_i U 4; 

12: rfe y — Ai fc Aj fc y; (Orthogonal projection) 

13: until ((||r fc || 2 > ||r fc _i|| 2 ) or (k > K)) 
14: k 4— k — 1; (Previous iteration count) 

Output: 

1: x G R w , satisfying xi^ = A^y and x^ fc = 0. 



introduce an input parameter 7 (0 < 7 < 1) that decides the 
cardinality of reduced set of potential atoms as V = L — \jL\ . 
In the iterations, the selection of atoms through orthogonal 
projection is performed in step [4] where V atoms are chosen 
from the set of L potential atoms. Then the chosen V atoms 



TABLE I 

Computational complexity comparison 



Algorithm 


Order of complexity 


Existing algorithms 


OMP 


0(K(MN + K' A + KM)) 


OLS 


0(K 2 (KN + MN)) 


SP 


0(K(MN + K 2 M)) 


Convex relaxation 


0{N 2 M) or 0(N 3 ) 


New algorithms 


POMP 


0(K(MN + K' A M)) 


LAOLS 


0(LK 2 (MN + K 2 + KM)) 


SOLS 


0((L - l-yL] )K 2 (MN + K 2 + KM)) 



are subjected through look ahead strategy such that a single 
atom is finally selected from the V atoms and added to the 
intermediate support-set. For 7 = 0, SOLS acts as LAOLS; 
on the other hand, for a 7 near to one that results in V = 1, 
SOLS acts as POMP. 

Each of the total matched filtering and the 
total orthogonal projection operations in SOLS 

« ELx [! + £'(* -(*-!))] = K+±L'{K*+K) = 
K + |(i - [jL\)(K 2 + K). Given the sparsity level K 
and the number of potential atoms L, the computational 
complexity decreases linearly with the increase in 7. For 
example, the choice of 7 = 0.5 leads to a reduction of 
computational complexity by nearly half compared to the 
LAOLS algorithm. 

V. Comparison of Computational Complexity 

In this section, we perform an analysis of computational 
complexity of the new algorithms and compare them with ex- 
isting algorithms. We mainly quantify the number of addition 
and multiplication operations. 

Let us first consider the POMP algorithm. For fc'th iteration, 
the algorithm requires one matched filter operation and two 
orthogonal projections. A matched filter operation requires 
2MN computations and each orthogonal projection requires 
0(k 2 M) computations. Overall the associated complexity is 
0(MN+K 2 M) where we assume the worst case scenario that 
k = K. Now, considering K iterations, the overall complexity 
of the POMP algorithm is 0(K(MN + K 2 M)). 

Next we consider the LAOLS algorithm. In this case, the 
intensive computation is due to the look ahead part where it is 
necessary to perform matched filter and orthogonal projection 
operations. We compute orthogonal projections recursively 
using the block-wise matrix inversion mechanism shown in 
Appendix [C] Using recursive computation and considering the 
worst case scenario that k = K, the necessary computation 
for each orthogonal projection is approximately 7K 2 +AKM. 
So, the total computation for each matched filter and each 
orthogonal projection is 0{MN + K 2 + KM). Now, consid- 
ering 7}L(K 2 + K) matched filter operations and orthogonal 
projections, the total complexity of LAOLS is 0(LK 2 (MN + 
K 2 + KM)). Following the similar arguments, the complexity 
of SOLS algorithm is 0((L - [jL\ )K 2 {MN + K 2 + KM)). 

We show the comparison of complexity between several 
algorithms in Table Q] Following [44] (pages 6-8 of [44]), we 
assume that complexity of convex relaxation methods is either 
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0(N 2 M) (using linear program) or 0(N 3 ) (using interior- 
point methods). The complexity of the existing IGS algorithms 
are evaluated in Appendix [A] and [B] Let us consider a 
quantitative analysis to show the complexity benefit of LAOLS 
over convex relaxation methods. From the literature, we note 
that most of the sparse reconstruction algorithms are successful 
when M = OiK log N). For our comparison, we assume that 
K varies linearly with N as K = tN (where < r < 1) and 
a typical setup where K = tN < M = KlogN < N. 
We also use L — K for LAOLS and find its computa- 
tional complexity as 0(K 3 MN); the term K 2 + KM is 
neglected as K 2 + KM <C MN. Therefore, the computational 
complexity of LAOLS is 0(K 3 MN) = 0(K 4 N log N) = 
0(t N 5 log N). We also assume that interior-point methods 
are used for convex relaxation algorithms and hence their com- 
plexity is 0(N 3 ). To achieve a better computational advantage 
for LAOLS than the convex relaxation algorithms, we must 
need C(r 4 iV 5 log N) < 0(N 3 ). This may be satisfied when 
T 4 /Y 5 log/Y < iV 3 , or r < (N 2 logiV) -0 - 25 . For example, 
when N = 500, we need r < 0.0283 and hence the allowed 
K < 0.0283 x 500 15. The above analysis is also valid for 
SOLS. 

To achieve complexity advantage for LAOLS/SOLS over 
convex relaxation algorithms, the above analysis shows the 
need of a lower K. The complexity increases heavily for a 
higher K and hence, it is natural to seek suboptimal solutions 
for a higher K. A possible suboptimal solution may be not 
to allow the full future evolution in look ahead strategy, but 
to allow some extent by restricting the number of future 
iterations. Considering a depth variable r\ that denotes the 
necessary allowable future iterations performed in look ahead 
strategy, a trade-off between performance and computational 
complexity may be achieved by proper choice of L and -q. 
An engineering intuition is that the increase of any of them 
keeping the other fixed leads to a saturation in performance. In 
the next section, we allow the look ahead strategy to check the 
full future evolution and show saturation in performance for 
increasing L. We also show running time comparison results 
between competing algorithms. 

VI. Experiments and Results 

We performed computer simulations in order to compare the 
performance of following CS reconstruction algorithms: OMP, 
OLS, SP, POMP, LAOLS, SOLS, and basis pursuit/basis pur- 
suit denoising/LASSO (BP/BPDN/LASSO) [10]-[11], [14]- 
[15]. Among these algorithms, BP/BPDN/LASSO is a convex 
relaxation algorithm (li norm minimization based) and the 
other algorithms are IGS methods. The BP/BPDN simulation 
code (matlab) was taken from the Zi-magic toolbox [18]. 
The LASSO simulation code (matlab) was taken from the 
SparseLab toolbox [19]. For SOLS, we used a fixed value of 
7 = 0.5. Experimentally we evaluated the effect of L on the 
trade-off between complexity and reconstruction performance 
for POMP/LAOLS/SOLS. We first discuss the reconstruction 
performance measures and experimental setups, and then 
discuss the BP/BPDN/LASSO algorithm briefly followed by 
reporting the performance results of all the algorithms for 
clean and noisy measurements. 



A. Performance measures and experimental setups 

We use two performance measures. For the first performance 
measure, we use signal-to-reconstruction-noise ratio (SRNR) 
defined as 

£{||x||I} 



SRNR 



£{\\*-niY 



(6) 



where x is the reconstructed signal vector. Note that our 
objective is to achieve a higher SRNR. 

Next we define another performance measure which pro- 
vides a direct measure of estimating the underlying support 
set. For a A"-sparse signal vector x, the support set was 
denoted as I with cardinality K. Let us denote the support 
set of reconstructed vector x as I. We assume that x is 
also a X-sparse signal vector, i.e. \I\ = K. To measure the 
support set estimation error, we consider to use the distortion 
d(I,l) = 1 — (\I n 2\/K^j [45]. Considering a large number 
of realizations (signal vectors), we can compute the average 
of d(I,I). We define the average support-cardinality error 
(ASCE) as follows 



ASCE = £ { d(l, 2)} = 1 - — 8 { \1 n i|} 



(7) 



Note that the ASCE has the range [0, 1] and our objective is 
to achieve a lower ASCE. Along-with SRNR, the ASCE is 
used as the second performance evaluation measure because 
the main objective of the IGS algorithms is to estimate the 
underlying support set. 

Now we discuss experimental setups. In a CS setup, all 
sparse signal vectors are expected to be exactly reconstructed 
if the number of measurements is more than a certain threshold 
[35]. However, the computational complexity to test this uni- 
form reconstruction ability is very high. Instead, for empirical 
testing, we can devise a strategy that can compute the perfor- 
mance measures for random measurement matrix ensemble. To 
measure the level of under-sampling, let us define the fraction 
of measurements (FoM) 

«-#■ 

Using a, steps of the testing strategy are listed as follows: 

1) For given values of the parameters K and N, choose a 
such that the number of measurements M is an integer. 

2) Randomly generate a sensing matrix A G R MxN 
where the components are drawn independently from a 
Gaussian source (i.e., a t j ~ Af (0, ij) and then scale 
the columns of A to unit-norm. 

3) Randomly generate a set of X-sparse data x where 
the support set I is chosen uniformly over the set 
{1,2,..., N}. Let we denote the size of data as S (i.e., 
the number of signal vectors x is 5"). The non-zero 
components of x are independently drawn by either of 
the following two methods. 

a) The non-zero components are drawn independently 
from a standard Gaussian source. This type of 
signal is referred to as Gaussian sparse signal. 

b) The non-zero components are set to ones. This type 
of signal is referred to as binary sparse signal. 
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Note that the Gaussian sparse signal is compressible in 
nature. That means, in the descending order, the sorted 
amplitudes of a Gaussian sparse signal vector's compo- 
nents decay fast with respect to the sorted indices. This 
decaying trend corroborates with several natural signals 
(for example, wavelet coefficients of an image). On the 
other hand, a binary sparse signal is not compressible 
in nature, but of special interest for comparative study, 
since it represents a particularly challenging case for 
OMP-type of reconstruction strategies [33], [35]. 

4) For each data, compute the measurement y = Ax + w 
and apply the CS reconstruction methods independently. 

5) Repeat steps 2-4 for a given times (let T times). Then 
evaluate the CS performance evaluation measures (by 
averaging over ST data). 

6) Repeat steps 1-5 for a new a. 

This test can be performed for any chosen K and N. 

Considering the measurement noise w ~ J\f (0, ct^Im), we 
define the signal-to-measurement-noise-ratio (SMNR) as 



the absence of a slack estimate e, we can solve an li norm 
penalized l 2 norm cost function in the following way: 



SMNR: 



gjMf> 
£{||w||I}' 



(9) 



where £{||w||2} = o^M. We report the experimental results 
at varying SMNRs. 

In the presence of a measurement noise, it is impossible 
to achieve perfect CS recovery. On the other hand, for the 
clean measurement case, perfect CS recovery of a sparse signal 
is possible if a exceeds a certain threshold. In the spirit of 
using CS for practical applications with a less number of 
measurements at clean and noisy conditions, we are mainly 
interested in a lower range of a where performances of the 
contesting algorithms can be fairly compared. 

B. BP/BPDN/LASSO 

In this subsection, we discuss BR BPDN and LASSO [10]. 
The BP/BPDN is a convex relaxation algorithm that looks 
for directly minimizing the li norm of the solution vector 
with respect to convex constraint functions. For the case of 
clean measurement (i.e. w = in fl}), we use BP where the 
following optimization problem is solved 

Pi : x = argmin ||x||i subject to y = Ax. (10) 



The above optimization problem can be recast as a linear 
program (LP) [18]. For the case of noisy measurements, we 
use BPDN [10]. Using a quadratic constraint, we solve for 



arg mm 

xGR" 



|x||i subject to ||y - Ax|| 2 < e, (11) 



where e > || w|| 2 is a parameter that depends on the noise prop- 
erties. The above problem can be recast as a second order cone 
program (SOCP) [18] and solved by an interior-point method. 
An interior-point method, such as the log-barrier method [44], 
is powerful and provides consistent performance. According 
to the suggestion of [12], we choose e 2 = a 2 w \ M + 2\/2M) . 
The ^i-magic toolbox [18] provides BP and BPDN codes. In 



oi 



argmin {A||x||i + ||y - Ax|| 2 } . (12) 
xeR N 



The above formulation is called LASSO where the regu- 
larization parameter A needs to be optimally chosen. The 
optimal choice of A can be performed by LARS (least angle 
regression) [15]. SparseLab toolbox [19] provides code for 
such an optimal LASSO. 

We note that BP/BPDN/LASSO does not provide a solution 
vector that is guaranteed to be if-sparse in any experimental 
condition. Therefore, to compute ASCE performance measure, 
we assume that the support set is constructed from the solution 
vector of BP/BPDN/LASSO by considering the atom-indices 
associated with K highest amplitude coefficients. Then a 
new reconstructed signal x is estimated using the orthogonal 
projection of y on the constructed support set of cardinality K 
(LS solution). We used the new reconstructed signal for our 
purpose of evaluating SRNR. In the literature [46], [47], such 
an oracle driven approach is shown to provide better SRNR 
performance for convex relaxation algorithms. Through ex- 
perimental evaluations, we also verified that the oracle driven 
approach improves the performance of BP/BPDN/LASSO 
significantly [47]. 

C. Experimental Results 

Using N = 500, K = 20, S = 100 and T = 100, we 
performed experiments. That means, we used 500-dimensional 
sparse signal vectors with sparsity level K = 20. Such a 4% 
sparsity level is chosen in accordance with real life scenarios, 
such as most of the energy of an image signal in the wavelet 
domain is concentrated within 2 — 4% coefficients [10]. We 
used 100 realizations of A (i.e., T = 100). For each realization 
of A, we used 100 signal vectors that are randomly generated 
(i.e., S = 100). Then, we incremented a from a lower limit 
to a higher limit in a small step-size (with the constraint that 
corresponding M is an integer for a value of a). Therefore, for 
each CS method at a chosen a, the performance is evaluated 
through averaging over 100 x 100 = 10000 realizations. 

Let us first observe the effect of increasing L for the 
POMP, LAOLS and SOLS algorithms. For the Gaussian sparse 
signal at clean measurement condition, we show the SRNR 
performance of the three algorithms in Fig. [T] We show the 
results for L = 2, 5 and 20, in the range of a from 0.1 
to 0.2. In Fig. Q] (c), we do not show the performance of 
SOLS for L = 2 (for L = 2, SOLS acts as POMP). From 
Fig. Q] we observe that the SRNR performance improves as L 
increases and the performance improvement by increasing to 
L = 5 from L — 2 is similar to the performance improvement 
by increasing to L = 20 from L = 5. This observation 
corroborates that the performance improvement saturates with 
increasing L. Considering the user defined choice that L < K, 
we choose to use L = 20 as it provides the best performance 
for all the three algorithms. For binary sparse signal, we also 
performed similar experiments and found the same trend in 
performance. In the later part of this subsection, we show 



9 




Fig. 1, Effect of increasing L on the POMP, LAOLS and SOLS algorithms for Gaussian sparse signal at the clean measurement condition. SRNR (in dB) 
versus fraction of measurements (a) is plotted: (a) POMP, (b) LAOLS, and (c) SOLS. 



the results for POMP, LAOLS and SOLS algorithms where 
L = 20 is used. 

Next, we compared between OMP, OLS, SP, BP/BPDN, 
LASSO, POMP, LAOLS and SOLS algorithms. We first show 
the SRNR and ASCE performance results for Gaussian sparse 
signal in Fig. (0. The results are shown for clean measurement 
condition and noisy measurement conditions with SMNR = 20 
and 10 dB. Fig. 12(a) and (b) compare the performance of all 
the algorithms for clean measurement condition. The SRNR 
results of Fig. [2] (a) is shown in Fig. [3] for a better view. We 
observe that the OLS and OMP provide similar performance. 
BP provides better SRNR than LASSO, but similar ASCE. 
The SP performs poorer than the OLS and OMP algorithms. 
For the proposed algorithms, we find that all the three new 
algorithms perform better than the existing algorithms. POMP 
provides a considerable improvement over OMP, and even 
performs better than the BP. At a = 0.2, the POMP provides 
more than 6 dB SRNR performance improvement compared 
to OLS and OMP algorithms. The LAOLS outperforms all the 
other competing algorithms by a significant margin and can 
be considered the best. At a = 0.2, the LAOLS provides 
nearly 20 dB SRNR performance improvement over OLS 
and OMP algorithms and 15 dB improvement over BP. The 
SOLS and POMP are found to be the second and third best, 
respectively. Now, let us observe the performance trends for 
noisy measurement conditions. Fig. [2] (c) and (d) show the 
comparative study at SMNR = 20 dB. In this case, we still find 
that the new algorithms are better than the existing algorithms. 
LAOLS and SOLS are found to show significant improvements 
than the others. A natural question is what happens if the 
noise power increases. At SMNR = 10 dB, Fig. [2] (e) and (f) 
show the results where we note that the BPDN and LASSO 
provide similar performance trends and they are better than 
the other IGS algorithms. For the Gaussian sparse signal, the 
results show that the IGS algorithms are non-robust at a higher 
measurement noise power. 

Now we show the comparative results for binary sparse 
signal in Fig. (0J. The results are shown for clean measurement 
condition and noisy measurement conditions with SMNR = 
20 and 10 dB. Fig. [4] (a) and (b) compare the performance of 
all the algorithms for clean measurement condition. We again 
observe that the OLS and OMP provide similar performance. 
However, in this case, the SP is found to perform better 
than the OLS and OMP algorithms. The convex relaxation 



algorithms (BP and LASSO) provide the best performance 
among all the algorithms. For the proposed algorithms, we find 
that POMP provides a considerable improvement over OMP. 
It can be said that LAOLS is the best performer among the 
IGS algorithms. Let us now compare the algorithms for noisy 
measurement conditions. At SMNR = 20 dB, Fig. |4] (c) and 
(d) display the results where the trends are similar to the clean 
measurement case results shown in Fig.|4](a) and (b). However, 
the performance trends take an interesting turn at a higher 
noise power level. At SMNR = 10 dB, Fig. |2](e) and (f) show 
the results where we observe that the BPDN/LASSO does not 
provide the best performance, but the LAOLS provides the best 
performance closely followed by the second best performance 
of SOLS. For the binary sparse signal, the LAOLS, SOLS and 
SP are found to perform better than the BPDN/LASSO at a 
higher measurement noise power. 

Comparing all the algorithms for two different input signals 
at varying noise power and number of measurements, we note 
that the proposed POMP, LAOLS and SOLS have a promise 
to provide good performance. They consistently outperformed 
OMP and OLS algorithms. A specific example is that, for 
the Gaussian sparse signal at clean measurement condition, 
LAOLS provides 20 dB SRNR improvement than the OMP 
and OLS at a = 0.2. 

Next we performed a running time comparison between the 
competing algorithms for the CS reconstruction of Gaussian 
sparse signal at SMNR = 20 dB. The running time means 
the average execution time. To compute the running time, 
we used matlab codes executing on a desktop computer and 
the averaging is performed using CS execution times for 
100 data (we used T = 10 ans S = 10). We performed 
the running time experiments for different (N, K) pairs. The 
chosen pairs are (500,10), (500,20), (500,30), (1000,20), 
(1000,30) and (1000,40). For each (N, K) pair, the choice 
of M = \K\ogN~\. The running time results are shown in 
Table ITT1 Following the discussion in section [V] we denote the 
allowable K by K* for which the LAOLS and SOLS might 
have a lower complexity than the BPDN complexity (BPDN 
uses interior-point method). From Table [TT] we find that the 
OMP has the minimum running time. The POMP and SP have 
similar order of running time. The running time of LAOLS and 
SOLS increases heavily with the increase of K. In this table, 
SOLS has the parameter 7 = 0.5. Depending on the parameter 
7, the SOLS running time can be adjusted. We performed the 



For Gaussian sparse signal 
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Fig. 3. Comparison of OMP, OLS, SP, BP, LASSO, POMP, LAOLS and SOLS algorithms for Gaussian sparse signal in clean measurement condition. 
SRNR (in dB) is plotted against the fraction of measurements (a). The same results are shown in Fig. [2] (a). 



TABLE II 
Running time comparison 



Algorithm 


N = 500 


N = 1000 


K* = 15 


K* = 20 


K = 10 K = 20 K = 30 


K = 20 K = 30 K = 40 


Existing algorithms 


OMP 
OLS 
SP 
BPDN 
LASSO 


0.0010 0.0023 0.0043 
0.3076 1.0100 2.6124 
0.0057 0.0193 0.0482 
0.4599 0.8140 1.1478 
0.0442 0.1325 0.2700 


0.0042 0.0108 0.0220 
1.0522 2.6973 6.0095 
0.0226 0.0541 0.1063 
2.6550 2.9915 3.5564 
0.2358 0.5379 1.1774 


New algorithms 


POMP 
LAOLS 
SOLS 


0.0098 0.0340 0.0945 
0.0752 0.6696 2.7874 
0.0337 0.2560 1.0172 


0.0407 0.1099 0.2465 
1.1794 6.4881 22.1076 
0.4395 2.3020 7.6874 



running time comparison for binary sparse signal also, but do 
not report the results as they are similar to the Gaussian sparse 
signal case. 

We end this section with interesting observations. We note 
that, depending on the signal statistics (Gaussian and binary), 
the algorithms show different performance trends. For ex- 
ample, LAOLS provides better performance for a Gaussian 
sparse signal in the clean and moderately noisy measurement 
conditions. On the other hand, BP/BPDN/LASSO provides 
better performance for a binary sparse signal in the same 
testing conditions. Another example is that the SP provides 
poor performance for the Gaussian sparse signal, but found to 
be good for the binary sparse signal. Therefore, a challenging 
research problem is to develop adaptive IGS algorithms which 
are robust to underlying signal statistics. 

Reproducible results: In the spirit of repro- 
ducible results, we provide necessary download- 



able matlab codes in the following website: 
https://sites.google.com/site/saikatchatt/softwares/. For CS 
reconstruction of Gaussian sparse signal at SMNR = 20 dB, 
the codes produce the SRNR and ASCE results shown in 
Fig. |2] (c) and (d), and the running time comparison results 
shown in Table [TTJ 

VII. Discussions, Conclusions, and Future Works 

For greedy algorithms where atoms are selected one-by- 
one serially over iterations, we show that new schemes can be 
developed for selecting an atom from a set of potential atoms. 
The potential atoms are found using a standard matched filter. 
The use of standard matched filter reduces the search space 
reasonably and allows us to use more sophisticated tools for 
atom selection. Among the developed tools, the look ahead 
strategy allows us to see the evaluation of an algorithm in 
future iterations. Overall, for selecting atoms, the new tools 
help to improve further on the state-of-art performance of 
matched filter. 

Using our experimental framework, we also note that the 
performance of greedy search algorithms depend on the input 
signals' statistics. Future investigations can be performed to 
develop new greedy search algorithms that are robust to the 
statistics of input signals. 
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Appendix 

A. Serial Atom Selection Based OMP and OLS 

In this section, we briefly describe two serial atom selection 
based IGS algorithms: OMP and OLS. First we summarize 
the main steps of the OMP in Algorithm Q (see [31], [32] and 
Algorithm 3 of [33]). Note that, in each iteration, the OMP per- 



Algorithm 7 : OMP for CS Recovery 
Input: 

l: A, y, K; 
Initialization: 

1: Iteration counter k <— 0; 

2: r <-y,lo<- 0; 
Iterations: 

i: repeat 

2: ki-k + 1; 

3: ik <— index of the highest amplitude of A t r / t_ 1 ; 
4: X k <- J fc _a U i k ; (Note: \I k \ = k) 

5: r-fc <— y — Ax k Aj fc y; (Orthogonal projection) 

6: until ((||r fc || 2 > ||r fc -i|| 2 ) or (k > K)) 
i: k <— k — 1; (Previous iteration) 

Output: 

1: x e M. N , satisfying Xx k = Aj fc y and x^ = 0. 



forms a matched filter operation and an orthogonal projection. 
Using block-wise matrix inversion, the orthogonal projection 
operations can be performed recursively. Considering the worst 
case scenario that k = K, the necessary computation for each 
orthogonal projection is approximately IK 2 + 4KM. So, the 
total computation for each matched filter and each orthogonal 
projection is 0{MN + K 2 + KM). Now, considering K 
iterations, the total complexity is 0(K(MN + K 2 + KM)). 

Next we summarize the OLS algorithm [31] in Algorithm[8] 
The overall structure of the OLS algorithm is similar to 



Algorithm 8 : OLS for CS Recovery 
Input: 

l: A, y, K; 
Initialization: 

1: Iteration counter k <— 0; 

2: r <-y,I <- 0; 
Iterations: 

i: repeat 

2: k <- k + 1; 

3: i k <- arg min i:X(u) =z k _ 1 ui \\y ~ A Z M A z (u) yh> 

(i (£ Zk-i an d orthogonal projections) 
4: X k <- 2fe_i U i k ; (Note: \I k \ = k) 

5: r fc <r- y - A Xk A^y; (Orthogonal projection) 

6: until ((||r fe || 2 > ||r fe _i|| 2 ) or (fc > K)) 
7: k <— k — 1; (Previous iteration) 

Output: 

1: x € M N , satisfying Xx k = A x k y an< i ~X-i k = 0. 



the OMP algorithm. The only and important difference is 



the way how the atoms are chosen through iterations. In 
each iteration, an atom is selected in such a way that its 
choice leads to the minimum residual norm for that iteration. 
This approach of atom selection is different from the OMP 
which uses a matched filter to select an atom, but does not 
guarantee to provide the minimum residual norm. However, 
OLS is computationally more intensive than the OMP. The 
OLS requires to perform approximately [N + (N — 1) + 
.,. + (N - K)\ = \{KN - K 2 + K) orthogonal projec- 
tion operations. As typically N ^> K, the computational 
complexity can be significantly high. However, in this case 
also the orthogonal projection operations can be performed 
recursively using block-wise matrix inversion. Considering the 
worst case scenario that k = K, the necessary computation 
for each orthogonal projection is approximately IK 2 +4KM. 
So, the overall complexity of OLS is 0{KN{K 2 + KM)) = 
0(K 2 (KN + MN)). 

B. Parallel Atom Selection Based SP 

In this section, we describe the parallel atom selection based 
IGS algorithm: SP. We summarize the main steps of the SP 
algorithm in Algorithm [9] (see Algorithm 1 of [35]). The SP 



Algorithm 9 : SP for CS Recovery 
Input: 

l: A, y, K; 
Initialization: 
1: Iteration counter k 0; 

2: Iq ^— indices of the K highest amplitudes of A*y; 
3: r <- y - A^A^y; 
Iterations: 
l: repeat 

2: k 4- k + 1; 

3: Z( p ) {indices of K highest amplitudes of A*rfc_i}; 

4: Z {u) 4- I fe _i U l {p) ; (K < \Z {u) | < 2K) 

5: xj (u) 4— Aj- } y; Xj ^ 0; (Orthogonal projection) 
6: Ik 4— {indices of the K highest amplitudes of x}; 
7: r fe <r- y - A Xk Aj fc y; (Orthogonal projection) 

8: until (||r fc || 2 > ||r fc _i|| 2 ) 

9: k <— k — 1; (Previous iteration) 

Output: 

1: x e TRL N , satisfying xj fc = A Xk y and x^ = 0. 



algorithm starts with an initial if -element support set Iq and 
an initial residual ro = y — Ai Aj o y. At the fc'th iteration 
stage, it forms the 'matched filter' A*rfc_i, identifies the K 
highest amplitude coordinates, forms a dummy support set 
I{u) = Ik-i U I(p), refines out A"-element support set Ik 
from I( u ), solves a LS problem with the selected indices in 
2fc, subtracts the LS fit and produces a new residual. Given 
the sparsity level K, the algorithm estimates a support set of 
cardinality K in each iteration and runs until the residual norm 
minimization condition is violated. Note that, unlike in the 
case of serial atom selection based OMP and OLS algorithms, 
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here the support set cardinality is not increased one-by-one 
through iterations. Rather, a K -element support set is refined 
through iterations by addition of potential new atoms and 
deletion of unnecessary atoms. An important point is to note 
how the i'T-element support set If. is chosen from the dummy 
support set T/ u -\ through using the orthogonal projection that 
invokes LS solution. The dummy support set Ir u ^ is formed 
through unionizing the previously estimated support set Ik-i 
with the set of K new atoms' indices. Then the observation y 
is orthogonally projected on the span of atoms that are indexed 
in Ity\ followed by picking up K indices corresponding to the 
highest amplitude coefficients of the solution vector. In the 
context of POMP/SOLS algorithm, we used similar strategy 
in Algorithm []} Algorithm [5] for improving the atom selection 
strategy where a single atom/multiple atoms is/are selected 
from a set of potential atoms in each iteration. Normally, the 
SP algorithm converges less than K iterations. Assuming the 
worst case scenario that the SP algorithm runs at most K 
iterations, it performs K matched filtering and 2K orthogonal 
projection operations. Note that, like OMP and OLS, the or- 
thogonal projections can not be performed recursively. Hence, 
the total complexity for each iteration is 0(MN + K 2 M). 
Therefore, considering K iterations, the overall complexity of 
SP algorithm is 0(K(MN + K 2 M)). 

C. Recursive Computation for Look Ahead Strategy 

In Algorithm [5] a computationally intensive part is per- 
forming the orthogonal projection operation in each itera- 
tion. For orthogonal projection, we need to perform pseudo- 



inverse of Aj fc e 

,t 



where Aj 



[A 2 



The 



At 

-J-k 



(14) 



Now, considering the fact that the multiplication of two ma- 
trices, C £ R mx ™ anc [ J) £ R nxfc , requires 2kmn operations, 
we can show that the total required computation for pseudo- 



inverse of Ax k £ 



pAf xfe 



is « 7k 2 + 4kM. 
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